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ABSTRACT 

Three  techniques  for  finding  the  eigenvalues  and  eigenfunctions 
are  investigated,   A  typical  problem  involves  a  linear  homogeneous 
differential  equation  with  variable  coefficients  of  the  form 

P(x)  y"(x)  +  p'(x)  y«(x)  +  a;2  M(x)  =  0  .  (l) 

The  functions  p(x)  and  M(x)  are  functions  which  are  positive,  or  have 
at  most  isolated  zeroes  on  the  fundamental  interval  (0,L);  tu  is  a 
parameter.   Appropriate  end  conditions  are  specified  so  that  the 
problem  is  self-adjoint. 

The  three  methods  are:  Newton's  method,  Stodola's  method,  and 
■khs  H-'vlei'0"^1™"?^ tz  method*  Tbr:  methods  ore  derived  and  a  T^TYn+.oT^ 
solution  by  each  method  is  included  in  the  paper.  A  second  problem 
involving  Bessel's  equation  of  order  zero  is  solved  using  each  method 
and  a  comparison  of  the  eigenvalues  and  eigenfunctions  is  made  with 
tabulated  values. 

The  results  indicate  that  Newton's  method  is  to  be  prefeired 
usually. 
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I.   INTRODUCTION 

Consider  the  linear  homogeneous  differential  equation  with 
variable  coefficients  of  the  form 

P(x)  y"(x)  +  P'(x)  y'(x)  +  u>2  M(x)  y(x)  =  0,    xc(0,L),      (1) 

where  P(x)  is  of  class  C  ,  M(x)  is  continuous,  and  "both  functions  are 
positive,  or  have  at  most  isolated  zeroes  on  the  fundamental  interval 
(0,L).   Also  given  are  end  conditions  of  the  form 

a1  y(x.,)  +  b1  y  (x^  =  0 

(2) 

a2  y(x2)  +  b2  y'(x2)  =  0 

The  function  y  and  the  parameter  uj  are  to  be  determined. 

An  equation  of  the  form  (l)  may  arise  as  a  result  of  solving  a 
partial  differential  equation  of  the  form 

o  t 

by  separation  of  variables  with  boundary  conditions  corresponding  to 
(2),   The  solution  of  (3)  is  assumed  to  be  of  the  form 

U(x,t)  =X(t)  T(t)  (4) 

and  standard  techniques  for  separation  of  variables  [8]  are  employed. 
Equation  (l)  may  also  occur  as  the  Euler  equation  [l]  of  an  isopermetric 

problem  of  the  form 


L 
J      P(x)   y'(x)      cbc  =  KEN  (5) 

0 


subject   to  the   constraint 

L 


2 


J       M(x)  y(xr  cbc  -  1    .  (6) 

0 

This  is  equivalent  to  solving  the  problem 

L  2     L 

P(x)  y'(x)  dx/J  M(x)  y(x)2dx  =  MIN.        (7) 
0  0 

In  this  case,  the  Euler  equation  is 

h  ^Fy«)  "  Fy  =  ^  [p(x)  y'(x)]  +  A(x)  y(x)  =  °        (8) 

where 

F  =  f °  -  co2  f1  (9) 


anu 


f°  =  p(x)  y(x)2 


J  =  M(x)  y(x)2. 


(10) 


Calculus  of  variations  treats  the  minimization  or  maximization 
of  functionals  such  as  integrals.   This  paper  considers  three  techniques 
for  obtaining  the  eigenvalues  and  eigenfunctions  of  Eq.  (l).  These 
techniques  are  Newton's  method,  Stodola's  method,  and  the  Rayleigh- 
Ritz  method  [3,  2,  5]. 

Newton1 s  method  may  be  applied  to  the  solution  of  either  the 
isopermetric  problem  or  the  equation  resulting  from  separation  of 
variables.  The  principal  problem  is  to  obtain  a  value  of  the  solution 
to  the  differential  equation  which  will  satisfy  the  end  conditions 
imposed  by  Eq.  (2). 


Consider  y  as  a  function  y(x,co)  on  those  extremals  [l]  from 
(x-l  t   y-i)'   Let  the  derivative  of  y  with  respect  to  co  be 

3  y  (x..u>) 
S(x,u>)  =—^3 .  (11) 

A  correctional  equation 

y(x2,"j)  -  y(x2) 

\ew  =  ^Old  "    sCx2,tu)     "  t12' 

is  derived  in  Chapter  II  and  the  sequence  of  steps  to  obtain  the 
solution  is  discussed  there. 

Now  consider  Eq.  (l)  written  in  the  self-adjoint  form 


_d_ 
dx 


P(x)  y'(x)   ==  co2  M(x)  y(x)  (13) 


subject  to  the  constraints  (2).   In  Stodola's  method,  consider  (13)  to 
have  the  form 

L  y  =  co2  y   .  (14) 

where  the  operator  L  is  defined 

"W^S^J'W)  (15) 

2 

For  the  eigenvalues  co.  and  associated  eigenf unction  Xi, 

L  Xi  =  u)f  Xi  .  (16) 

Taking  into  account  the  end  conditions,  define  an  inverse  operator, 
L   ,  then 

L"1  Xi  =  \   Xi  (17) 

CO. 

1 


Since  differential  equation  (13)  and  the  boundary  conditions  (2) 
are  self-adjoint,  there  is  an  infinite  set  of  orthogonal  eigenfunctions, 
Xi,  and  an  arbitrary  function,  Y~,  satisfying  the  end  conditions  (2), 
can  be  expanded  in  a  series  of  eigenfunctions 

YQ  =  a1  X1  +  a2  X2  +  ...  (18) 

Repeated  application  of  the  inverse  operator  (17)  to  Y~  emphasizes 
the  coefficient  of  X-,,  and  decreases  the  relative  size  of  the  other 
coefficients.  The  resulting  function  Y  is  the  corresponding  eigen- 
function.  Further  details  of  Stodola's  method  are  discussed  in 
Chapter  III. 

The  Eayleigh-Ritz  method  is  a  procedure  for  obtaining  approximate 
solutions  of  variational  problems.   The  procedure  consists  of  assuming 
that  the  desired  extremal,  y(x)  can  be  approximated  by  a  linear 
combination  of  n  suitably  chosen  functions 

y(x)  «  C1  ^(x)  +  C2  $2(x)  +  ...  +  CM  $M(x)  ;  (19) 

The  C .  are  to  be  determined  to  effect  the  desired  minimum.  Usually 
the  $  i(x)  are  chosen  to  meet  the  boundary  conditions  for  any  choice 
of  C,  An  approximation  of  y(x)  may  also  be  made  by  spline  functions 
of  the  form 
y(x)  ~  C1  +  C2  x  +  C3  x2  +  C4  [(x  -  0^)]  +  ...  +  C^  [(x  -  o^)2 J 

where  the  C.  are  to  be  determined  and  the  o/.  are  in  (0,L).  The 
1  x        v 

function  (x  -  a  )       =  (x  -  a-  )2  j_f  x  >  a     and  is  zero  for  x  <  o^  . 


The  C's  are  not  independent  but  must  "be  chosen  to  satisfy  the  end 
conditions.   In  both  the  polynomial  and  spline  function  approximations, 
determination  of  the  C.  is  accomplished  by  solving  an  algebraic 
eigenvalue  problem  [4]. 

In  the  following  chapters,  Newton's,  Stodola's  and  the  Rayleigh- 
Ritz  mentod  are  dismissed  in  more  detail  and  computer  solutions  of 
(l)  by  the  three  methods  are  presented. 


II.      NEWTON'S  i-ETHOD 

The  problem  of  minimizing  the  numerator  of  Eq.  (1#7)  with  the 

denominator  constrained  to  equal  one  is  again  being  considered.  The 

problem  may  be  expressed  in  the  form 

L  L  n 

I  =  f  P(x)  y'(x)2  dx  =    f  dx  =  MIN  (1) 

0  0 


and 


L  L 

J  =  J  M(x)  y(x)  dx      f  dx  =  1  (2) 

0  0 

Assume  the  desired  solution  y(x)  =  y  (x)  has  been  found.   Then 

(5) 
a2  y(x2)  +  b2  y  (x2)  =  0 


and 


1  -x-        •  -x-' 

J  f  (x,y,y')  ^  =  1  ,   y  =  y    ,   y1  =  y  (4) 

0 

A  curve  satisfying  conditions  (5)  and  (4)  is  called  admissable. 


A.   DERIVATION  OF  NECESSARY  CONDITION 

The  first  necessary  condition  is  derived  for  the  case  where 
b.  =  0  =  bp  .   Consider  replacing  y  (x)  by  y  (x)  +  el]  (x),  where 
Tl(x)  is  an  arbitrary  function  vanishing  at  0  and  L  and  having 
piecewise  continuous  first  derivatives.   It  is  necessary  that 

g-£($-J^ax)n*Max--o  (5) 


10 


or 

L 


g  =  J     MQ(x)   Tlf(x)   dx  =  0  (6) 


0 
for  all  T|  such  that 
L 


-I     (fl<    -  I     f1v  **)   V(x)    dx  =  0  (7) 


dJ 

cle  "  JQ   -y'        JQ  ' 


or 

L 


dJ 
de 


~  =  J     M  (x)   T\   (x)    dx  =  0  (8) 

0 


In  Eq.    (6)    and   (8),   M. (x)  is  defined  as 

x 

M.(x)    =   (f1,    -    P  f1  dx)            *              ,          „, 

ix  '      v  y'     Jn  y      y  -  y     >    y    =  y* 


(9) 


Let 

L 
g  Ax)   =  M  (x)   -  VL        =  M  (x)   -  J     M  (x)    dx/L   .  (10) 

1  1  av         x  0 

The  above  necessary  condition  may  be  reexpressed.      It   is  necessary 
that 

§  =  /"g  (x)  Vto   dx=0  (11) 

y         0 

for  all  T|  such  that 
t, 


dJ 

de 

0 


=  J     gl(x)   7]    (x)    dx  =  0  (12) 


since 


L  L  L 

J  Miav  1|  (x)  dx  =  Miav  J  T|'(x)  dx  =  Miav  Tl(x)  J  =  0       (13) 

Now  choose 

Tl'(x)  =  gQ(x)  +  c  Sl(x)  (14) 


11 


t  .  . 

and  choose  c  so  that  Tj  (x)  is  admissable: 

dT    L 

^  =  J  S^B0   +  c&1)  dx  =  0  .  (15) 

If  g  (x)  /-   0,  Eq.  (16)  can  be  solved  for  c  giving 

L  L      2 

c  =  -  J  g0(x)  g  (x)  dx  /  J  gQ(x)   dx  (16) 

0  0 


dl 

de 


=  J  s0(^)   (s0(x)  +  o  g.,(x))  dx  (17) 


Now  consider 

L 

0 
then 

L 


+  °  de"  =  J*  (g0^X^  +  C  S1(x))(g0(x)  +  c  g^x))  dx         (18) 


dl  .    dJ 
de 


Then 


0 

(g^(x)  +  c  g„  (x)  )  dx  ^  0  . 

0  N  °         :    ' 


f>0  (19) 

unless 

gQ(x)  +  c  gl(x)  =  0  .  (20) 

Hence  this  is  a  necessary  condition,  known  as  the  First  Necessary 
Condition,  or  the  integrated  form  of  the  Euler  equation.  Writing 
Eq.  (20)  as 

f°.  +  c  f1,  -  f  (f°  +  c  f1)  dx  =  (lUx)   +  c  ML(i))  (21) 

yi      yi   Jq  v  y      y^      V  0W      1   yav 

0      1 

it  is  seen  that  f  and  f  enter  only  on  the  combination 

P  =  f°  +  c  f 1  .  (22) 
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The  integral  form  of  Euler's  equation  [l] 

x 
F  ,  -  f  F  dx  =  M   =  CONSTAM1  (23) 

yl     UQ   y         aV  V  Ji 

is  obtained  from  (20)  by  setting 

Mav  =  M0   +  c  M1    '  (24) 

av    uav      'av 

Now  let 

c  =  -  w2  (25) 

then  the  Euler  Equation 

Py.y.       y'W      +      Fy,y     ^  to      +     ^y  x        *(*)       "     ^y      =      0  (26) 

becomes 

P(x)  y"(x)  +  Pf(xj   y*(x)   =  -  uu2  M(x)   y(x)    .  .(2?) 

It  may  be  noted  that  this  differential  equation  is  self-adjoint. 

B.   COMPUTATIONAL  ROUTINE 

It  will  generally  be  necessary  to  generate  the  curve  by  forward 
numerical  integration.  There  are  two  parameters 

Cj  =  c  =  -  co2  (28) 

and 

c2  =  y'(0)  (29) 

which  are  needed  to  determine  a  curve  which  is  an  approximation  to 
the  solution;  an  estimate  of  y?(o)  is  needed  in  order  to  start 
integration. 


13 


In  order  to  apply  Newton's  method  it  is  necessary  to  evaluate 
derivatives  with  respect  to  these  parameters.  Let 

u .  |^M  (30) 

and 

Then  differentiating  Eq,  (27)  with  respect  to  c.  gives 

P(x)  u"  +  P  (x)  u'  -  c  M(x)  u  =  M(x)  y(x)  (32) 


Denote  the  left  side  of  (32)  as  Lu,  defining  the  second-order 
linear  differential  operator  L.   Then  the  equations  for  u  may  be 
written 

L  u  =  M(x)  y(x)  (33) 

and 

u(0)  =  0 

.  (34) 
u'(0)  =  0 

The  corresponding  equations  for  v  are 

P(x)  v  +  Pf(x)  v'  -  c1  M(x)  v  =  0  (35) 


and 


v(0)  =  0 


(36) 


v  (0)  =  1 


Ik 


or 

L  v  =  0 

v(0)=0  (37) 

v'(o)  =  1  , 

The  correctional  equation  for  J  is  obtained  from  the  approximate 
relations 

M  ■  H7 A  ci +  H; A  c2  (38) 

=  \   -  J  J 

by  J  is  meant  the  computed  value.   The  correctional  equation  for  y 
is  obtained  from  the  linear  approximation 

A  y(L)  =  u(L)  A  Cl  +  v(l)  A  c2  (39) 

=  -  y(L) 

where  y(l)  is  the  computed  value. 

However  this  problem  has  a  peculiarity  that  makes  it  somewhat 
simpler.  Because  both  integrals  are  quadratic  and  homogeneous  in  the 
pair  y(x)  and  y  (x) ,  the  principal  problem  is  to  find  c .  so  that 
y(L)  =  0.  The  constant  cp  is  then  obtained  by  substituting  the 


solution  y  into  J: 


c2  =  1/  *f~J  (40) 


Consider  then  y  as  a  function  of  a  single  parameter 

y  =  y(x,  cj  (41) 
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The  correctional  equation  reduces  to 


A  y  -  u(L)  A  Cj  (42) 

=  y(L)/u(L) 


or 


A  oj2  =  -y(L)/u(L)  (43) 

2 

In  the  computational  routine,  an  initial  estimate  of  ou  is  made. 

Then  y(x)  and  u  are  obtained  by  integrating  Eq,  (33) »   The  correction 

2  2 

for  u>  is  then  obtained  from  Eq.  (43)  •   If  a  good  guess  for  to  is 

made,  this  method  converges  rapidly.   If  a  normalized  value  of  y(x)  is 
desired  then  a  new  solution  may  be  obtained  from 


L  „    X2 


i 
2 


y(x)new  =  y(x)/(J  M(x)  y(x)<  dxj  (44) 

To  tell  if  the  lowest  eigenfunction  has  been  obtained,  a  check 
is  made  to  see  if  the  eigenfunction  obtained  has  a  zero  in  (0,L)« 
If  it  has  no  zero  in  (0,L),  the. lowest  eigenfunction  has  been 
obtained. 


16 


III.   STODOLA'S  METHOD 
Consider  Eq.  (1.1) 


_d_ 
dx 


(p(x)  y'(x))  =  -w2  M(x)  y(x),  (l) 

subject  to  the  constraints 

a1  y(0)  +  b1  y'(0)  =  0 


(2) 


a2  y(L)  +  b2  y'(L)  =  0 


Stodola's  method  of  finding  the  lowest  eigenvalue  and  eigenfunction 

will  now  he  developed.   The  system  of  Eqs,  (l)  and  (2)  is  self-adjoint, 

2 


* 


i  =  1,  2,  , ..  and  associated  eigenfunctions  X, 
Consider  Eq,  (l)  as  having  the  form 


L  y(x)  =  oj  y(x)  .  (5) 


in  which  the  linear  differential  operator  L  is  defined  by 

and  consider  any  particular  eigenfunction  Xi,   Then 

L  Xi  =  cu2  Xi 
1 

L2  Xi  =  vfi     Xi  (5) 


Tm  v.    2m  ^. 
L  Xi  =  w.   Xi 

1 
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_1 

The  inverse  operator,  L   ,  is  now  defined,  taking  into  account 

the  end  conditions  (2),  Consider 


_d_ 
dx 


(p(x)  Xi*  J  =   -  cu?  M(x)  Xi 


(6) 


and  integrate  (6)  to  obtain 

P(x)  X±'  =  -  wi  (j  M(x)  Xi  dx  -  c) 


(7) 


the  constant  c  is  to  "be  determined  later.  The  eigenf -unction,  X. , 
may  be  obtained  by  dividing  Eq,  (7)  by  P(x)  and  integrating  again 


X.  =  -  U3 

1       1  L 


.  X 

I. 


M(x)  Xi  dx 

0 

— KZ 


X 


dx 


J0  *&T 


dx  -  d 


(8) 


The  constants   of  integration,    c  and  d,   may  be   obtained  from  the 
ftones-fcrair-te?    in   Rn      ( ?)      and    F,ns      (7)    ^^^    ''p^         ^nimno.  +Vipqp> 


equations   gives 


c  = 


J     M(x)   X±  dx 
^2  I  P(x) dx  + 


b     J     M(x)   X     dx 

0 
P(L) 


dx 


b 


2 


a2  \ 


(9) 


and 


r»  ax  <l  d.     j 

La2  J0   HxT +  pTlj  -  i~m\ 


d  =  -  b^  d/a1  P(0) 


(10) 


if  a1  ^  0.   If  a,  =  0  or  P(o)  =  0  then  c  =  0  and  a  solution  for  d 
is  obtained. 

The  set  of  operations  on  the  right  side  in  Eqs.  (7)  and  (8)  define 
the  inverse  operator,  L~  ,  with  c  and  d  given  by  Eqs.  (9)  and  (10), 
These  operations  give 
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X.  =  uj2  L"1  X.  (11) 


or 


l"1  x.  =  x./^2  (12) 

1     1!     1  v  ' 

Stodola's  method  makes  use  of  the  fact  that  an  arbitrary  piecewise- 
continuous  function  which  satisfies  the  end  conditions  (2)  can  be 
expanded  in  a  series  on  the  eigenf unctions.  The  first  estimate  Yo 
of  the  first  eigenf unction  may  be  chosen  as  follows.  Let  Yo  be  the 

first  estimate  of  X.:   choose  it  so  that  it  has  no  zeroes  inside  the 

1 

interval  (0,L).  Then  Yo  may  be  considered  to  be  the  following  form 

Yo  =  a     X±  +  a2  X2  +  ...  (13) 

where  a.,  is  not  equal  to  zero;  For  convenience  Yo  is 'normalized;  the 
L^  norm  of  Yo,  |  |Yo| |,  is  defined  as 

MAX  |Yo(x) I  =  | |Yo| |  .  (14) 

C  <  X  <  L 

and  I  I Yo  j   is  set  equal  to  one.   Now  consider 

L~1Yo  m  -1  X  +  -|  X  +  ...  (15) 

1       2 


and 


m1       m2 


uu 


2m         uk  2m 


1  c- 
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From  Eq.  (16)  it  can  be  seen  that  the  relative  size  of  the  components 
Xp,  X.,,  . ..  is  decreasing  by  a  factor  ( — )   ,  ( — J  ,  , .,  respectively. 

-1      ^ 

After  a  few  applications  of  the  operator,  L   ,  to  Eq.  (13),  the 
leading  term  will  dominate.   It  is  convenient  to  normalize  after 
each  iteration;  let 

Z  =  L"1  Y  ,  (17) 

m       m-1  v  ' ' 

1 
where  Ym  is  the  m  th  approximation  to  X,    and 


Ym  =  Zm/    MZmM    •  (18) 

When  this  is  compared  with  Eq.  (12),  where  X.  on  the  left  hand  side 
corresponds  to  Y  ,  and  X.  on  the  right  hand  side  corresponds  to  Y  , 


it  is  seen  that  tu  is  approximately 


2 


»1  -1/llzJI  •  (19) 

The  iteration  is   stopped  when 

HYmW  -VlWH    <e   '  (2C) 

where  e  is  seme  specified  small  number.   The  resulting  function  Y 
is  an  approximation  to  X. . 

In  the  computational  routine  an  initial  approximation,  Y  to 
X.  is  made  which  satisfies  the  constraints  (2)  and  Y  is  normalized. 
The  next  approximation  to  X„  is  obtained  by  forward  integration  of  Eqs, 
(7)  and  (8).   The  cont rants  of  integration,  c  and  d  arc  determined 
by  Eqs.  (9)  and  (10).   The  value  of  u>  is  approximately  1/||Z  ||. 
The  new  approximation,  Z  ,  is  normalized,  Eq.  (18),  th-  procedure  is 
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repeated.   When  the  norm  of  the  difference  between  successive  iterates, 

Y  . .  Y  ,  is  less  than  some  prescribed  value  the  routine  is  stopped. 
m-1 '  m7  ** 

The  value  of  cu  from  Eq,  (19)  is  taken  as  the  lowest  eigenvalue  and 

Y  is  the  corresponding  eigenfunction.   The  numerical  solution  of 
Bessel's  equation  of  order  zero  gave  a  good  approximation  to  the 
eigenvalue  but  a  smaller  value  of  the  eigenfunction. 
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IV.   RAYLEIGH-RITZ  METHOD 

The  Rayleigh-Ritz  method  has  long  "been  used  to  obtain  an  approximate 
solution  to  eigenvalue  and  eigenf unction  problems.  This  is  still 
a  useful  method  if  a  high  speed  computer  is  not  available. 

Assume  that  m  functions  § .  (x)  are  given  which  satisfy  the  end 

conditions.   These  are  often  chosen  as  polynomials  or  trigonometric 

functions.   The  solution  then  is  approximated  by  a  linear  combination 

of  these 

m 
y  =  S  c  S  (x)  ,  (1) 

i=1 

This  is  substituted  into 

?   J*  i   ?  k         ? 

uT  --    P(x)  y  (x)   dx  /    M(x)  y(x)   dx  =  MIN         (2) 

0  0 

and  the  c's  are  chosen  to  minimize  (2).   This  is  equivalent  to  the 
minimization  of  a  quadratic  form 

m    m 

S    £  a. .  c.  c .  , 
i=1   i=1   ^  1  J  ' 


where 


L       , ,\  . 


■j   ~0 


subject  to  a  constraint  that 


a.   =  I     P(x)  §.(x)  6  (x)  dx  ,  (3) 


m   m 

E   S   b.  .  c.  c.  =  1  ,  (4) 

i=1  >1   1J  x  3 
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where 


L 


b   =  J  M(x)  3  (x)  I  (x)  dx  .  (5) 

-J    0       i     j 

Equation  (2)  may  then  be  written 


co2  =  C  A  C  =  fflN  (6) 

C  B  C 

where 

1=  a    ,   i  =  1,2,,.,,  m  ,   j  =  1 ,2, , . , ,m,  (7) 


B  =  b±..   ,   i  =  1,2,,,,,  m  ,   j  =  1,2,.,,,m,  (8) 

and 

_  m 

C  =  (c1f  c?,  *».  c  )  (9) 


The  problem  may  be  reduced  as  follows.   Let  the  eigenvalues  of 

2 

B  be  L  ,  i  =  1,2,,,,ra  and  the  associated  normalized  eigenvectors 

u .  ,   Then 
l 

?BT  =  I  (10) 

where 

T  =  (u^  u2,  ..,,  u  )  BIAG  (1/X1f  1A2,  ...,  1A  )         (11) 

and  the  transformation 

C  =  T  D,  (12) 

where 

_  T 

I)  =  (d1?  d2,  ..,,  dm)  (13) 
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gives 


-jji  —       _rp  _m  _____ 
D     D  =  D     T     B  T  D 

-T--  (14) 

=  C     B  C 


=  1 

The  desired  minimum  value  of 


(?  A  C  =  1?  T^  A  T  D  =  ICL'IT  (15) 


subject  to  (14)  is  the  minimum  eigenvalue  of 

m—    h   m 


(16) 


and  the  corresponding  vector  D  is  the  associated  normalized  eigen- 
vector of  E.   Since  the  elements  of  T  and  D  are  known,  C  may  "be 
evaluated  from  (4.12), 

The  minimization  technique  using  the  spline  function  approximation 

y(x)  «  C1  +  C2  X  +  C3  X2  +  C4  [(x  -  ^)2]  +...+  CK+3  C(x  -  cy2] 

(17) 

is  similar  to  the  polynomial  approximation  method.   In  this  case 
however,  tne  functions  $.(x)  do  not  satisfy  the  end  conditions. 
Hence  the  constants  above  are  not  independent  but  one  of  them  is 
dependent  on  the  others  for  the  end  conditions  to  be  satisfied.   It 
is  necessary  to  have  an  initial  subprogram  to  eliminate  this  constant. 
Then  the  problem  is  reduced  to  the  same  form  as  before 


m-1  m-1  ,   , 

£   2  a.  .  c.  c.  «  TCLN  (18) 

i-1  .  «   x3  *  J 

11  0=1 


2k 


and 


m-1  m-1 
Z       £  b.  .  c.  c.  =  1  (19) 

i=1  j=1  -J  ±    J 


and  solved  as  before. 
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V.   APPLICATION  OF  RAYLEIGH-RITZ , 
STODOLA  AMD  MEV/TON  METHODS 

In  Chapters  II,  III  and  IV,  three  methods  were  developed  for 
finding  the  lowest  eigenvalue  and  eigenfunction  for  Eq,  (1.1).  These 
methods  were  used  to  obtain  numerical  solutions  of  the  following 
two  problems:  Bessel's  equation  of  order  zero 

X  y"(x)  +  y'(x)  +  w2  x  y(x)  =  0  ,   xe(0,2.405),  (1) 

subject  to  the  constraints 

y(0)  =  1  ,   y(2.i+05)  =  0  ,    y'(0)  =  0  ;  (2) 

and  an   equation  of  the   form 

(1+x)   y"(x)   +  y'(x)   +  uj2  y(x)    =  0    ,        xe(0,l)    ,  (?) 

subject  to  the  constraints 

y(o)  =  y(i)  =  o  .  (4) 

In  this  chapter  the  computation  and  numerical  results  are  discussed. 
In  the  Rayleigh-Ritz  polynomial  approximation,  an  approximation 
for  the  extremal,  y(x) ,  satisfying  Eq»  (3)  was  made  by  choosing 
suitable  functions  3>  (x)  and  §?(x)  to  satisfy  the  constraints  in 
Eq.  (4)»   The  approximation  for  y(x)  is 

y(>0  =  o1  ^(x)  +  c2  $2(x)  (5) 

2 

=  C.  x  (1-x)  +  c2  x   (1-x) 
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The  values  of  c. ,  c _,  and  the  lowest  eigenvalues  were  obtained  by 
the  techniques  described  in  Chapter  II.   Computer  Output  1.  shows 
the  values  of  c  and  c?,  the  lowest  eigenvalue  and  the  values  of 
y(x)  at  tenths  of  an  interval.   The  numerical  solution  of  Eq.  (3) 
using  the  spline  function  approximation 

y(x)  «  c.,  x  +  c2  x2  +  c^  [(x-1/3)2]  +  c4  [(x-2/3)2]  (6) 

is  shown  in  Computer  Output  2.  and  Computer  Program  2.   This  solution 
followed  a  format  similar  to  the  polynomial  approximation  solution. 
However,  a  transformation  was  required  to  eliminate  the  dependent 
coefficient.   In  this  case  the  coefficient  c  was  taken  as  dependent 
on  c., ,  c„  and  c.,.   The  required  matrix  products,  transpositions  and 
eigenvalue  solutions  were  evaluated  by  using  subroutine  OMPRD,  GMTRA, 
and  EIGEN  respectively,  [4]. 

A  numerical  solution  of  Eq.  (3)  by  Stodola's  method  using  as  a 
first  estimate 

y(x)  =  4x  (1-x)  (7) 

to  satisfy  the  constraints  in  Eq.  (4)  yielded  the  values  in  Computer 
Output  3.   The  solution  of  Eq,  (3)  by  Newton's  method,  with  an  initial 
estimate  of 

uj2  =  2.0  (8) 

is  shown  in  Computer  Output  4»   ln  these  two  methods  the  fundamental 
interval  (0,1)  was  sub-divided  into  one  hundred  equal  sub-intervals 
to  perform  the  second  order  Runga-Kutta  integration  routine.   If  the 
error  between  the  computed  value  and  the  boundary  value  at  X  equal 
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one  was  less  than  one  ten  thousandth,  the  program  was  terminated. 
The  programs  for  Stodola's  and  Newton's  method  are  shown  in  Computer 
Program  3*  and  Computer  Program  4» 

The  lowest  eigenvalues  obtained  by  the  above  methods  are  reasonably 
close  in  value.  The  maximum  difference  in  values  was  between  Newton's 
method 


au  =  3-78634 


and  Rayleigh-Ritz  polynomial  approximation, 


w  -  3-79313  . 


(9) 


(to) 


Newton's  method  converged  in  three  iterations  and  Stodola's  method 
required  seven  iterations. 

The  Ra"'j  ei ~\i~ P. it 7. ,    Stodola  and  Newton  methods  were  vsed  to  find  a 
numerical  solution  to  vessel's  equation  of  order  zero,  Eq,  (l),  A 
tabulated  solution,  [7l»  of  Eq.  (1)  gives  the  smallest  eigenvalue 
as  one.  The  result  obtained  by  the  Rayleigh-Ritz  method  was 


u>  =  1.01435 


(11) 


when  the  function  was  approximated  by  quadratic  splines 


y(x)  =  c„  +  c2  X  +  c_ 


(x  -  0.8)* J  +  c  [(x  -  1.6) 


(12) 


The  smallest  eigenvalue  of  Eq,  (1)  obtained  by  using  Stodola's 

method  with  an  initial  estimate  of 

.2 


y(x)  -  ((1-x)/2,405J 


(15) 
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IS 

uj  =  1.00000  (14) 

In  Newton's  method  an  initial  estimate  of 

«J  =  1.4  (15) 

yielded  a  value 

uj  -  1.0000  (16) 

for  the  lowest  eigenvalue. 

A  comparison  of  the  eigenf unction,  y(x) ,  obtained  by  the  three 
methods  was  made  with  tabulated  values  for  y(x)  of  Eq.  (1).  The 
values  obtained  by  Newton's  method  at  tenths  of  an  interval  were 
wd thin  five  ten-thousandths,  by  Rayleigh-Ritz  were  within  five 
thousandths  and  by  Stodola's  method  were  within  two  one-hundredths . 
The  best  approximation  of  the  lowest  eigenfunction  was  obtained  by 

Newton's  method, 

2 

The  solution  of  ou  obtained  by  the  Rayleigh-Ritz  method  is 

greater  than  the  theoretical  value  as  given  in  Reference  5»  The 

numerical  solution  of  Bessel's  equation  of  order  zero  supported  this 

2 

conclusion.  The  value  obtained  for  w  was  approximately  one-thousandth 

larger  than  the  theoretical  value. 

The  results  of  the  three  methods  are  shown  in  Computer  Output 
5.,  6, ,  and  7»  and  the  programs  are  shown  in  Computer  Programs  5»» 
6. ,  and  7» 
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VI.   CONCLUSION 

The  three  methods  each  yield  satisfactory  values  for  the  eigen- 
values and  the  eigenvectors.  Generally,  Newton's  method  seems  to  be 
most  satisfactory.   It  is  for  the  most  part  straightforward  to  program. 
The  only  significant  errors  that  are  not  apparent  are  there  due  to  the 
integration  routine.   It  is  very  flexible,  being  applicable  to  many 
different  types  of  problems.   It  converges  rapidly  if  a  good  initial 
estimate  is  made,  and  computing  times  seem  to  be  generally  small. 
Proof  of  the  exact  convergence  of  Newton's  method  may  be  found  in 
Reference  3* 

Stodola's  method  also  is  quite  satisfactory.   There  seems  to  be 
■/■us   Qiiestior  o^  coziv^rfrencs   Th"-  initial  =cucss  can  be  csuite  nocr 
and  it  will  still  converge.  More  iterations  than  for  Newton's 
method  seems  to  be  required  to  get  the  same  accuracy,  and  there  is  no 
simple  way  to  get  an  estimate  of  the  error. 

The  Rayleigh-Ritz  method  seems  to  have  little  to  recommend  it 
if  a  good  computer  is  available.   It  is  more  difficult  to  program. 
It  has  one  advantage  that  no  iteration  is  required.  Use  of  spline 
functions  increases  the  tedium  of  programming  so  that  they  are  not 
worthwhile. 

The  methods  of  Stodola  and  Rayleigh-Ritz  were  used  before  the 
advent  of  large  scale  computers,  though  application  was  rather  limited 
compared  to  problems  that  can  be  treated  now. 
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In  case  larger  eigenvalues  and  associated  eigenfunctions  are 
needed,  Uewton's  method  yields  these  directly  if  the  initial  estimate 
of  uj  is  changed.  The  number  of  zeros  of  the  eigenfunction  obtained 
show  which  eigenvalue  and  eigenfunction  ha.s  been  obtained.  The 
eigenf unctions  are  automatically  orthogonal  with  a  weight  factor 
M(x)  since  the  system  is  self-adjoint. 

If  Stodcla's  method  is  used  to  obtain  higher  eigenvalues  and 
eigenf unctions  it  is  necessary  to  adjoin  a  further  condition  that 
each  new  eigenfunction  be  orthogonal  to  all  of  the  preceeding.  This 
increases  the  tedium  of  programming  and  decreases  the  accuracy, 
since  preceding  values  are  only  approximate. 

In  the  Rayleigh-Ritz  method,  an  estimate  of  other  eigenvalues 
and  eigenf unctions  is  obtained  from  the  other  eigenvectors  of  the 
matrix  S  of  chapter  IV.  The  eigenvalues  are  too  high  unless  the 
eigenfunction  can  be  expressed  as  a  finite  sum  of  the  functions 
$.(x). 
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COMPUTER  OUTPUT  1.      (Polynomial) 

THE    SMALLEST    EIGENVALUE  ?o 798129 

THE    VALUES    OF    CI     AND    C? 

Cl=  7oC142cl 

C2=  -3ol385?0 

THE    VALUES    OF    Y(X)     AT    CME-TENTH    TNTEPVAL 


X 

Y 

OoC 

0, 

n 

0o 1COCCO 

0- 

603^39 

Oo 200000 

lo 

021852 

0 ,300000 

1 

275272 

•5,4000?  " 

1 

382  >  3~ 

%  500or° 

1 

3612?6 

^o 59P999 

1= 

231482 

Oo  699999 

1: 

011638 

°o799999 

0- 

72C556 

no  8999°Q 

0 

3  77r6  8 

C0  909999 

0 

OOOC    4 
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COMPUTER  OUTPUT   2.       (Spline) 

THE    SMALLEST    EIGENVALUE  3„ 79C108 

VALUES    OF    C1*C2,C3,AND    C4 

Cl=  6o358921 

C2=  -7ol97122 

C3=  0o565650 

C4=  5o281205 

VALUES    OF    Y{X)     AT    ONE-TENTH    INTERVAL 


X 

Y 

GoO 

OoO 

OolCOOCO 

Oo 5  63920 

Oo  200000 

0o983899 

Co  300000 

10259933 

0o4000C0 

lo394543 

Oo500COO 

lo 3  95  892 

Do  599999 

lo  264613 

Co  699999 

lo006572 

Co  799999 

0o698055 

Oo 89999  9 

0o362533 

0o999999 

Oo 00000  3 
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COMPUTER  OUTPUT   3.      (Stodola) 

NI=    7cOMEGA    =  3o78~7142 

THE    MAXIMUM    DIFFERENCE    BETWEEN    YB    AND    YBN    IS  OoOD0047 

PE,KE,AND    ROOT    (PE/KE)     ARE  7>172187  r0 498484  3o 793123 

X  YBN 

0  o  r  0  o  c 

OolOOCCO  Do  387621 

0o2C00CC  n0 6942 88 

Cp3C0C0C  Oo 897866 

0»4C0.CO0  0o992039 

no5C0CCn  Oo 982  3  33 

O06OGCC-O  0o 383094 

Co  700000  Oo 713473 

A-800rfr.  f>0  49  5322 

0o90CCC0  0o25f609 

loooooro      -O0OOCD14 


3h 


COMPUTER  OUTPUT  k.      (Newton) 

THIS    IS  PATH    NOo     3o     OMEGA=                  3,78634 

X  Y 

^o^  OoO 

OolOCOO?  Qor>93143 

Co200C0C  Co  166838 

0o300C0r  r0?i576  8 

Cc4000CC  Oo 238397 

Co50000C  C3 236090 

%6C0C?C  ^o 212242 

ro7C0CDC  ro 17149  8 

O08OOCOC  Do  119087 

0o90000C  0o 06^299 

lo 000000  O0OOOO8I 
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COMPUTER  OUTPUT   5.      (Spline) 

THE    SMALLEST    EIGENVALUE  lo0001C3 

VALUES    OF    C1,C2,C3,AND    C4 

Cl=  -U132401 

C2  =  0o272184 

C3  =  -0ol25127 

C4  =  -Oo 184540 

VALUES    OF    Y(X)     AT    ONE-TENTH    INTERVAL 


X 

Y 

AoO 

-lol32401 

Oo 240500 

-lo 116657 

Oo 481000 

-lo 069427 

Oo 721500 

-Co990712 

Oo 962000 

-Oo 88379^ 

lo 20  249  9 

-Oo 7  59094 

lo 443000 

-0o617383 

lo  683499 

-Oo 459947 

1,924000 

-Oo 302293 

2o 164499 

-0ol48986 

2 o 405000 

-0o0000C9 
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COMPUTER  OUTPUT  6.      (Stodola) 

NI=    7c    CMEGA=  lcCOCOOC 

THE    MAXIMUM    DIFFERENCE    BETWEEN    YB    AND    YBN    IS  OoOnOC21 

PE, KE, AND    ROOT     (PE/KE)     ARE  C. 755341  0o761183  0.996155 

X  YBN 

OoC  le OOO CCO 

0.2405C0  Go 9 740 24 

0.481000  0.931927 

0o72150C  0« 8637  87 

0*962000  Co  772 530 

lo 202499  Co  662059 

1.4430C0  0.537061 

lo 683499  Co  402789 

U 924000  Co  2  64  808 

2.164499  Co  128741 

2.405000  OoO 


37 


COMPUTER  OUTPUT  7.   (Newton) 


THIS  IS  PATH  NO.  4 
X  Y 

0.0  1.000000 

0.2*10500  O.985594 

0.1+81000  0.91+2999 

0.721500  0.874051 

0.962000  0.781713 

1.202499  0.669932 

1.1+43000  0.543451 

1. 683499  0.407585 

1.924000  0.267966 

2.164499  0.130282 

2.405000  0.000010 


OMEGA  = 


0.99998 
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COMPUTER    PROGRAM    lo 

RAYLEIGH-RITZ    METHOD    APPLIED    TO    EQN    (lol) 
(l+XJ^'MXI+YMXI  +  l  W**2)*Y(  X)=t. 


DIMENSION    PH(2f?.r  1  )  , PHP(2,  101  )  ,A(2,2t  ,3(2,2),  P<101)  , 
lFA(!'l),FBM.HtBL(2f2)  ,AMB(2  ,2  )  , C  (2  )  ,  AX  (3  )  ,  BX  <  3  )  , 
2RLB<2),  BXD(2)  ,  T(2,?),D(?  ) ,TT(2,2)  ,  ALPHA (2, 2) , AP(2,2) , 
3S(2t2) «RLBD(2) tY(lCl) fZ(lOl) 

REAL    LEV 

X  =  ro' 

DX=0.-C! 

DP  1  1=1  I '  1 
C  DEFINITION  0^  ELEMENTS  OF  MATRIX  B 

P(I)=1    C+X 

PH<1, I )=X*( loP-X) 

PH(2,  I )=(X**2) *(loO-X) 
C    DEFINITION    OF    ELEMENTS     OF    MATRIX    A 

PHPdiI  )  =  lo"~2 


T-*X-3oO*(X*,*2) 


AND    B(I,K) 


PHP(2, I )=2 

x=x+(>0'  1 

1    CONTINUE 
COMPUTE    A( I , K) 
DO    2    1=1,2 
DO    ?    K=l,2 
FAINT    =n 
FBINT    =r 
FAI1  )  =  P<1 )*PHP ( 
FB(1)=PH(I, 
90       a    j  =  ?  ,  1  ■  i 

FA( J)=P( J)*PHP( I , J)*PHP(K, J) 
FB( J)=PH(I, J)*PH(K, J) 


!  1  ' 
'  H  ( K  * 


*PHP(K, 1 


FAINT  =  FAINT    +(FA(J)+    F  A(  J-l  )  ) ''DV? 

FBINT=FBINT+    ( FB ( J ) +  FB( J-l )  H  DX/" 

CONTINUE 

A(  I,K)=FAINT 

B(  I  TK) -FBINT 


2^2 
201 


,2,H 

EVAL    O1 


3  CGNTINUF 
2  CONTINUE 
INITIALI7E    BX 

NC=1 

DO    ?f>  1    J  =  l,2 

DO    2."  2    1=1,2 

BX(NC)  =  B(  I» J) 

NC=NC+] 

IP     (IoGEoJ)    GO 
CO  NT  INl'E 
CONTINUE 
CALL    El  GEN     (BX, 
RLB( I)     IS    1/SQPT    OF 
KK=1 

DO    8    11=1,2 
BXD(  I  I  )  =  p.X!KK) 
KK=KK+I I+i 
CONTINUE 
DO    84     11=1,2 
DO    9    JJ=1,2 
B    (II, JJ  )  =  "\ 
CONTINUE 
CONTINUE 
DO    6    1=1,2 

RL3D( I )=SQRT(BXD( I )  ) 
RLB(I )=lo0/PLBD(I  ) 
B  (I, I)=PLB( !) 


FOR  THE  INPUT  TO  E1GEN  ROUTINE 


TO  2"1 


8 


9 
84 


39 


6  CGNTTNnr 
C    MULT       S(I,J)^:B    (  I,  J)=T<  I,  J  ) 

CALL    GMPRD     (S,B     ,T,2,2  ,2) 
C    TRANSPOSE    T=TT 

CALL    GMTPA    (T,TT,2,2) 
C    MULT       A*T=AP 

CALL    GMPRD     (     A , T t AP , 2, 2 , 2 ) 
C    MULT       TT*AD=ALPHA 

C&LL    GMPRD    (Tj  ,AP,  i\LPHA,  2,2,2  ) 
C    FIND    E    VAL    AND    E    VEC    OF    ALPHA 

NC  =  1 

DO    210    J =1,2 

DO    211    1=1, 
AX(NC) 

NONC+3 

IF    (  IoGEoJ)    GO    TO    210 
211     CONTINUE 
210    CONTINUE 

CALL    EIG^N    (AX  , S , 2 , ^ ) 

LEV=SQRT( AX (3) ) 

DO    7    1=1,2 

D(  I)=    S(  1,2) 

7  CONTINUF 

C    MULT       T*D    TO    GET    VALUES    nF    C 

CALL    GMPRD     ( T , D, C , 2 , 2 , 1 ) 
C    EVALUATE    Y    AND    7 

Z(l  )=0 

DO    1?    I=2,lr! 

Z(  I)  =  Z(  I-l)+ro     1 

12  CONTINUE 

DO    1?    I  =1  ,1     1 
Y(I)=C(1)*PH(1,I )+C(2)*PH(2,I) 

13  CONTINUE 

WRITE     (6,1     )     LEV,C(1)  ,C(2)  ,(Z(I  )  ,Y(  I)  ,1  =  1,1'  1,1**  5 
1C     FORMAT     (  '1' ////////, f*V  ,15X, » THE    SMALLEST    EIGENVALUE 
1F15c6//  ,  «C  ,  15X,  'THE    VMJJlS    OF    CI    AND    C2V/,«     •  ,15X, 
2,Cl=,,F15o6//,,C » ,  15X,»C2=' ,F15o6//, »C» ,15X, 
3»THE    VALUES    OF     Y(X)     AT    ONE-TENTH    INTERVAL1/', 
4,^« ,18X, »X» ,15X, »Y»//( '  « , ?F15o6//) ) 

wnr   (6,9°fi) 

°98     FORMAT     Cl'^X) 
STOP 
END 


i+O 


COMPUTER    PROGRAM    2o 

RAYLEIGH-RITZ    METHOD    APPLIED    TO    EQN    (loll 
(1  +  X)*Y«  '  (X>+Y«  (X  )+(W**2)*Y(X )=r 


DIMENSION    A  (3,3) 
1T(3,3), ALPHA(3 


B(3, 


3) , AX (6) 

AP(4,4)    ,TT(3 


2BXDC' )  ,  RLBD(3)  ,  AO(  4  ,4  )  .  RP  ( t- ,  u  ) 


,BX(6) 
t3) »S(3 

•  CM  (4 


,R(3,3),RLB(3) 
,3),D(3),C(3), 
, 3) ,CMT(3,4) ,BCM( 


34,3)  ,  ACM  (A,  3)  ,  P  (in  )  ,  PH(A,  1C  1  )  , PHP(^,lr  1  )  ,  FAQr  1  ) 
4,FB(101 ),Y(1     1) ,Z(1C1) 
REAL    LfV 

X  =  r'o 

DO    1    1=1,1'   T 
DEFINITION    OF    ELEMENTS     OF    MA'RIX    B 

P  ( I)  =  1 o  + X 

PH(1, I)=X 

PH(2, I  ) =X**2 

PH(3,T  )=((  X-l  o/3o  +  A.8S(X-lo  /3o  )  )/2o)*^2 

°H(4, I )= ( (X-2o/3 ,+ABS( X-2o /3o )  )/2o )**2 
DEFINITION  OP  ELEMENTS  OF  MATRIX  A 

PHP( i, 1 )=1 

PHP(? , I )=2  *X 

PHP  (3, 1  )=X-lo/3o+ABS(X-lo/3o ) 

PHP (4. I )=X-2o /3o+ABS( X-2o/3o ) 
92    X  =  X-K  o'  1 

1  CONTINUE 
DX=0,  01 

COMPUTE    AD(I,K.1     AND    BP(I,K) 
DO    2    I =]  ,4 

DO    3    K  =  l,4 
FAINT=T o 
FBINT='  o 

FA( 1 )  =  P ( 1  5  -  PHP ( I f 1 ) *PHP( K f 1 J 
FB(1 )=PH( 1,1 )*PH(K,1 ) 
DO    4    J -2.1'   ] 

FA( J)=P( J)  <PHP( !,  J )*PHP(K, J) 
FB( J)=PHf I ,  J)* ph(K, J) 
FAINT=FAINT+(FA(  J  )  +F  A(  J-l )  )  *DX/2 
FEINT  =FBINT  +  (F3(J)+PB(J-.l))-DX/2 
4    CONTTNUE 

APU,K)=FAINT 
BP(I,K)=FBINT 
3    CONTINUE 

2  CONTINUE 
INITIALIZE    Z* 

DO  51  KA=1 ,4 
DO  52  K3=1T3 
CM(KA ,KB)=C ,T 

52  CONTINUE 
51    CONTINUE 

EVALUATION  OF  THE  TRANSFORMATION  MATRIX  CM 
DO  53  KA=1 ,3 
CM(KA,KA)=lo' 

53  CONTINUE 
CM( 4, 1 1 --9, 
CM (4, 2 )=-9 
CM (4, 3) =-4 

CAL'  GMTRA(CM,CMT,4,3) 
MAKE  TRANSFORMATION  CM  (TR ANSPOSED)#B*CM 

CALL  GMPRDtB'5,  CM, BCM ,4,4,3 ) 

CALL  GMPRD(CMT, BCM,B,3,4,3) 
MAKE  TRANSFORMATION  CM ( TRANSPOSED) *A*CM 

C ALL    GMP RD ( A P , C  M , AC M ,4 , 4 , 3 ) 

CALL    GMPRD(CMT , ACM 

NC=1 

initialize;   bx  for   the 


,  A  ,  3 ,  4 ,  3  ) 

INPUT    TQ    EI  GEN    ROUTINE 


DO    2'  1    J=l,3 


111 


COMPUTER    PPOGRAM    3o 

STODOLA'S     METHOD    APPLIEP    TO    EQN(lol) 
(1+X)*Y»  «  (X  )+Y»  (X )+(W**2 )*Y(X)=r 


YB    IS    THE    GUESSED    VALUE    OE    Y    AND    YEN    IS    COMPUTED    VALUE    OFY 
DIMENSION    YB(1*1),SS(1    1 ) , P I ( 1C 1 ) , YBN( 1 r 1 )  ,  E R ( 1C 1 ) , 
1  X  { 1  0 1  ) 
REAL    KF 

1  CONTINUE 

S=^„ 

PI( l)=~o 

NI  =  1 

DX=^- ^1 

BFRO^='   o 

DO    15    1  =  1,1/1 

X(i )=0r i*(i-i) 

15  CONTINUE 

DO    2     1  =  1, lr  1 
SS(  D=r 

YE( I )=40*X( I )* ( 1-X( I) ) 
NORMALIZE    YB 

IF     (ABS(YE( I) ) dLEoBFROG)    GO    TO    2 
BFROG=ABS(YB( I ) ) 

2  CONTINUE 

DO    16    1=1,1"! 

YB( I )=YB< I  //PFROG 

16  CONTINUE 
?    C  t  •  N  F  T  N I  i  F 

DO    ',    1  =  1,1 
SO    =S 

S=S+(  YP  (  I) +vp(  i  +  i  )  )*DX/2c 

SS(I+1)=SS(T )  +  (SO/(?.+X(I ))+S/(l+X(I+l) ) )*DX/2o 
PI  (  1  +  3.  )  =  (!    /(I  +X(  I  )  )+lo/(l+X(  1+1)  )  )*DX/2,+PI(  I) 
^    CONTINUE 
NOW    GFT    NORM    OP    THE    NEW    VB 
PFROO    o 
DO    5     1=1,1 '1 

YBN(I)=PT(I)*SS(1^1)/PI(101)-SS(U 
IE     (ABS(YBM( I) ) oLE, BFROG)     GO    TO    5 
BFROG=ABS(YBN(  T )  ) 

5  CONTINUE 

OMB  =  SQRT(l-  /BFROG) 

ERR='o 

DO    6    1=1,1    : 

YBN(  I)=VF'N(  I  )/BF  ROG 

ER( T)=vb( I3-YBNM  I  ) 

IF    (ABS (ER( I ) loLT,  ERR)     GO    TO    6 

ERR=ABS(ER(  I  5 ) 

6  CONTINUE 
KE=C  o 

PE  =  r  o 

7  CONTINUE 

CHECK    TO    SEF    HOW    STODOLA    COMPARES    WITH    RAYLEIGH    QUOTIENT 
DO    R     1  =  1,1.  ' 

KE  =  K.p+(  YBN(  T  )  **?+YBN(  1+1  )**2)*DX/2 
PE=PE+(1+X( 1+1 ) )*< (Y3N(I+1 )-vBN( I ) >/DX)**2*DX 

8  CONTINUE 
OMB2=SORT{ PE/KE) 

IF     (NI0  LT07)     GO    TO    12 
Q    WRITF(6, K  )    NI  ,OMB, ERR ,PE,KE,0ME2, ( X(T  )  < VBN(I ) , 

1 1=1 ,  mi,  io) 

r-    FORMAT     (  •!»///////»  "*'  ,15X,»NI=»  ,I2,»o0l     GA    =',F12o6, 


h3 


12 


11 


270 

999 
998 


1/15X,«THF    MAXIMUM 

2F12o6,/15X, »PE,KE, 

3/2~X,'X'  ,lf  X,  'YEN' 


CONTI: 
DO  11 


OIFFEPENCE  EETWEFN  YB 
AND  ROOT  (PE/KE)  ARE' 
,//.  (12X,2F12o6// )  ) 


AND  YEN  IS' 
,3F12. 6,/ 


1  =  1,1 


YB( I) =YBN(I ) 

CONTINUE 

NI=NI+1 

IF  (NIr  LTo2?) 
GO  TO  999 
IF  (ERPoGTofo 
GO  TO  999 
WRITE  (6,998) 


GO  TO  27' 


)  )  GO  TO  3 


FORMAT 

STOP 

END 


(  '  1 '  ,  2  X  ) 


kk 


COMPUTER  PROGRAM  4, 

NEWTON'S  METHOD  APPLIED  TO  EQN  (Id) 
(l  +  X)vY»  « (X)+Y« (X)+(W**2 )*Y(X)=C 


C 

c 


GET 
USE 

100 


20< 


DIMENSION  X{ 101  )  ,Y( 101) ,YP< 101)  ,Y0( 101)  ,YP0(101) 
DO  1  1=1,101 

yd  )  =  o'"-i*(  i-i) 

CONTINUE 
Y(l)=Oo 

YP(l)=lo 
YFN=YP( 1  ) 
DX=OoCl 

NI  =  1 
YO  ( 1  )  =C 
YP0(1)=0 

YPON---YPOU  ) 

MODE  AND  OMEGA  BY  NEWTONS  METHODo 
RUNGA  KUTTA  2ND  ORDER  INTEGRATION 
0M=2o0 
DO  2  1=1 tlC« 
J=-l 

XN=X( 1+1 ) 
YN  =  Y( I ) 
YON=YO<  I  ) 

YFN  =  YP( I  )-(0M**2)*( Y( I )/( 1+X(  I )  )+YN/( 1+XN)  )*DX/2 
1-(YP (I )/(l+X(I ) )+Y°N/( 1+XN) )*DX/2 
YN=Y(  T  )+(YP(  T  )  +YPN  >*DX/2 
YPHN-YPO  in-ivcnm/iuvini  +YPDN  /  M  +XN  )  ) 


FIRST    GUESS 
ROUTINE 


OM 


:oV 


( 


( 1  +  X  ( I 


»     *    JU> 


ON/  {  1+XN)  )*DX/2-2~UM; 


v  t 


300 


400 


2+YfM/(  1  +  XN)  )*DX/2 
YON=YO( I )+(YPUl I )+Y£>0N)*DX/2 
IF    (JoGToO)    GO    TO    4< 
J  =  l 

GO    TO    20C 
YP( 1+1 )=YPN 
Y(I+1)=YN 
Y0(I+1 ) =  YON 
YPO( I+I ) =YPON 
CONTINUE 

IF     (NIo  LTo3)    GO    TO    11 
WRITE     (6,1. 
FORMAT     ('l1 


DX/2-(UM**-2) 
i  !  /  (  1  +  X  (  I  )  ) 


10 

"l«o     OMEGA 
11    CONTINUE 

NI=NI+1 
IF    ADMISSABLE     QUIT 
IF     (Y(Kl)oLTo  ( 
CORRECT    OMEGA 

0M=0M-Y( 101 )/Y0( 101) 
IF    NIoGTo2  5    QUIT 

IF     (NIo  LTo2?)     GO    TO    100 
500    CONTINUE 
STOP 
END 


)     NI,0M,(X(I),Y(I),I=1,1C1,1M 
,////////, «<j«  ,15X,«TH1S     IS    PATH 


•,F12o5,//,l3X,'X«,ilX,,Y',//( 


NOo ' , 12, 
:X,2F12o6//) ) 


0001)  )    GO    TO    50- 


h5 


COMPUTER  PROGRAM  5o 
RAYLEIGH-RITZ  METHOD  APPLIED  TO 
BESSEL'S  EQUATIUN  OF  ORDER  ZERO 


C  DEFI 


C  DEFI 


92 
1 

C  COMP 


6^2 


601 

600 
C  INIT 


52 
51 
C  FVAL 


53 


C  MAKE 


DIMEN 

T(3,3 

BXD(3 

BCM(4 

FB(10 

REAL 

X=0o 

DO  1 

P(  I  )  = 

EMC  I  ) 

NITIO 

PHC  1, 

PHC2, 

PHC3, 

PHC4, 

NITIO 

PHPC1 

PHPC2 

PHP(3 

x=x+c 

CONTI 
DX=(.o 
UTE  A 
DO  2 

DO  3 
FBINT 
FB(l) 
DO  4 
F  B  (  J  ) 
FBINT 
CONTI 
BPC  I, 
CONTI 
CONTI 
DO  6f. 
DO  60 
FAINT 
FAC1  ) 
DO  60 
FA(  J) 
FAINT 
CONTI 
API  It 
A(I,K 
CONTI 
CONTI 
IALIZ 
DO  51 
DO  52 
CMCKA 
CONTI 
CONTI 
UATIO 
DO  53 
CMCKA 
CONTI 
CMC1, 
CMC1, 
CMClt 
CALL 

TRAN 
CALL 
CALL 
NC  =  1 


SIGN  AC3, 3) ,8( 3,3) ,AXC 6)   ,BX(6)   ,  R ( 3  ,3 )  ,  RLBC 3 ) , 

) , ALPHA (3,3), AP(4,M,TT(3,3),S(3,3),P(3),CC3), 

)  , PLBDC3)  ,AQC4,4) , BPC 4, 4)  , CMC  4, 3)  ,CMT(3,4) , 

,3) , ACM (4, 3), PC 101) ,PH( 4,101) ,PHP(4, 101) , FAC101), 

1 ) ,EM(1l1 ) ,YC1C1) ,ZC1C1) 

LEV 

1=1,101 

X 

=  x 

N  OF  ELEMENTS  OF  MATRIX  B 
I  1  =  1. 

I )=(CX-o8+ABS(X-o8) )/2o )**2 
I)  =  C( X-1.6+ABSCX-1.6) )/2. )**2 

N  OF  ELEMENTS  OF  MATRIX  A 
, I >=2o*X 

,  I  )=X-DS+ABSC  X-o 8) 
.  I  )-X-lo6  +  ABS( X-1.6) 
.4  240  5 
NUE 
02405 
PC  I ,K) 
1  =  1,4 
K  =  l,4 


AND  BPC I,K) 


i  \ 

)   )*DX/2 


=  EMC1 )-PH(  1,1  )*dh(K,1 ) 

J=2?10J 

=  E  M  C  J )  -'  pH  <  T  -  • )  -"  pH  '  *  , 

=FBTNT+(FB( J )+FB( J-l 

NUE 

K)=FBINT 

NUE 

NUE 

0  1=1,3 

1  K=l,3 

=P(1)-PHP( I ,1 )*PHP(K, 1) 

2  J=2,l  I 

=  P( J)*PHP< I  ,J)*PHPCK, J) 
=  F  A  I  NT  +  {  F  A  (  J  )  +  F  A C  J- 1 )  )  * DX /  2 
NUE 

K)=FAINT 
)=AP< I,K) 
NUE 
NUE 
E  CM 

KA=1 ,4 

KB=1,3 
,KB)=Oo<" 
NUE 
NUE 
N  OF  THE 

KA=1,3 
+1 ,KA)=1.C 
NUE 

l)=-5o784 
2)=-2o576 
3)=-0o 648 
GMTRACCM,CMT,4,3) 
SFORMATION  CMC  TRANSPOSED) *B*CM 
GMPRDC BP,  CM,BCM,4 ,4,5 ) 
GMPRDCCMT,  BC'i ,  B  ,  3  ,  4,  3  ) 


RANSFORMATION  MATRIX  CM 


46 


C    INI 


202 
2C1 


C    RLB( 


8 


9 
10 


6 
C    MUL 


, R 1 3 , 0 ) 

o )     GO    TO    8f 6 
OF    Eo     VALo     OF 


C  TRA 
C  MUL 
C  MUL 
C    FIN 


TIALIZE    BX    FOR 

DO    201    J=l  ,3 

00    20  2    1=1,3 

BX(NC)=B(I,J) 

NC=NC+1 

IF    ( IoGFoJ )    GO    TO    201 

CONTINUE 

CONTINUE 

CALL    EIGEN     ( BX 

IF(BX(6)oLToC« 

I)     IS    1/SQRT 

KK=1 

DO    8    11=1,3 

BXD( I  I )=3X(KK) 

KK=KK+II+1 

CONTINUE 

DO    1C     11=1,3 

DO    9    JJ=1,3 

B    (II,JJ)=Go 

CONTINUE 

CONTINUE 

DO    6    I    =1,3 

RLBD( I )=SQRT(BXD(  I  )  ) 

RLB( I )  =  loC/RLBD(I  ) 

B     (I , I)=RLB( I) 

CONTINUE 
TIDLY  R( I ,J)"B  ( 

CALL  GMPRD(R,B 
NSPOSE  T=TT 

CALL  GMTRA(T,TT 
TIPLY  A*T=AQ 

CALL  GMPRD( A,T,A0,3,3,3) 
TIPLY  TT*AQ=ALPHA 

CALL  GMPRD(TT,A0,ALPHA,3 


THE  INPUT  TO  EIGEN  ROUTINE 


I, J)=T(  I, J) 
,T,3,3,3) 

,3,3) 


3,3) 


EoVALo  AND  EoVECc  OF  ALPHA 


DO  210  J =1,3 

DO  211  1=1,3 

AX(NC)=ALPHA(I , J) 

NC=NC+1 

IF(  IoGFoJ)     GO    TO    21C? 
211    CONTINUE 
210    CONTINUE 

CALL  EIGEN  (AX    ,S 
C  LEV  IS  THE  LOWEST  EI  GIN 

LEV=SQRT( A<(6)  ) 

DO  7  1=1  ,3 

D(I)=S(  I  ,3) 
7  CONTINUE 
C  MULTIPLY  T*D  TO  GET  VALUES  OF 

CALL  GMPRD(T,D,C  ,3,3,1) 

C0NE=-5o784*C( 1 ) -2o 5  76*C( 2 )-C 

Z(l)=Oo 


,3,0) 

VALUE 


C 


64G*C( 3) 


DO  64  1  =  2,  lv  1 


64 


1  ) 


Z( I)=0oC2405*(  I 

CONTINUE 

DO  6  5  1=1, 1 f  1 

Y(  I)=C0NE-PH(1, I )+C(l)*PH(2, I)+C(2)*PH(3,I)+C(3)* 


1PH(4,I ) 
65  CONTINUE 
WRITE  (6,1 
11=1, 101,1 
109  FORMAT 
1F15o6// 
2,'C1=  « 
3F1506// 
4« VALUES 
5'U' , 19X 
806  STOP 
END 


■Q)LEV,CONE,CU)  ,C(2),C(3),(Z(I),Y(I), 


( ' V  iff!!///,  »0«  ,15 
,  "J1  ,15X,  'VALUES  OF 
,F1536//,  "V  ,15X,  «C2=«  ,F15o6//, 


SMALLEST  EIGENVALUE1, 


C1,C2,C3,AND    C4'//, 


15X 


OF 


5X,'C4=     ',F1536//,» 


15X 


,15X,'C2 


Y(X)     AT     ONE-TENTH    INTERVAL'// , 


'X' , 14X, 'Y'//( « 


',2F15o6//J ) 


U7 


COMPUTER  PROGRAM  *D 

STCDOLA'S  METHno  APPLIED  TO 
BESSEL'S  EQUATION  OF  ORDER  ZERO 


C  YB  IS  THF  GUESSED  VALUE  OF  Y  AMD  YBN  IS  COMPUTED  VALUF  OFY 
DIMENSION  YBU  "l  )  ,  SS(  101  )  ,  PI  (101  )  ,vbn(]  ri) ,ER<1"1), 
1X( 101 ) 

REAL  KE 

1  CONTINUE 

S  =  ^o 

NI=1 

DX=<\  024r  5 

BFROG='  o 

DO  15  I=l,ir 1 

X< 1 )=C, 024^5* ( 1-1) 

15  CONTINUE 

DO  2  1=1,1' 1 

SS(I  )=' 

YB(I)=(l-X(I)/2o4C5)**2 

C  NORMALIZE  YB 

IP  ( ;,BS(  VB(  I  )  1  oLEoEPROG)  GO  TO  2 
BFRDG= ABS(ve( I ) ) 

2  CONTINUE 

DO  16  1  =  1,1' '  1 

YB< I)=YB( I J/BFROG 

16  CONTINUE 

3  CONTINUl 
S  S  < 1  *;  =  o 
SS(2)=YB(2KDX/2 

S=  'o 

S=S+(YB(1)*X<1  )+X(2  )*YB(2) )*DX/2 
DO  4  1=2,1  C* 

SO  =s 

S=S+(VP ( I )*X< I )+y< 1+] )*YB< i+i ) }*DX/2 

SSU  +  l)=SSU)  +  (  SO/X(I)+S/X(I  +  l  )  )*DX/2 

^  CONTINUE 
C  NOW  GET  NORM  OF  THF  NEW  YB 
BFROG='  o 
DO  5  1  =  1,  in 
VBNU  >=SS(1"  j  )-SS(  I  ) 
IF  (ABM  YBN(  I  )  )  oLEo  BFRO^)  GO  TO  5 
BFROG=ABS(YBN(  I ) ) 

5  continme 

OMB  =  SQRT(l.:  /BFPOG) 
EPR  =  r  o 

DO  6  i=i,m 

YBN(  I  >  =  YBN( I )/BFROG 

ER( I)=YB(  I)-YBM(  I  ) 

IF     (ABS(FR(  I)  )0LT0EPt?)     GO    TO    6 

ERR=ABS <ER(  I  )  ) 

6  CONTINUE 
KE=<so 

P  E  =  C1  o 

7  CONTINUE 

C    CHECK    TO    SEE    HPW    STOOOLA    COMPARES    WITH    RAYLFIGH    QUOTIENT 
DO    8     1  =  1,  K' 

KE=KE+(  X(  I  )^  YPN  (  I  )*  *2+X<  1  +  1  )*YBN(  I  +  i  )**2  )*DX/2 
PE=PE+X(I )*< ( YBN(I4  1)-YBN( I ) ) /DX)**2*DX 

8  CONTINUF 
0M32=SQRT(PE/KE) 

IF    (NIoLT07)    30    tp    1J» 

9  WRITF(6,1C]    NI ,OMB,ERR,PE,KE,OMB2, (X( I) ,v      '(I), 
11=1,101 ,10) 


)+8 


K-    FORMAT     ( '1'///////,  s     '  ,15X,*NI  =  ' ,12,  •<,    OMEGA=«,F12o6, 
1/15X,»THE    MAXIMUM    DIFFERENCE    BETWEEN    YB    AND    YBN    IS' 
2F12o6,/15X,f  <>E,KEt  AND    ROOT     (PE/KE) 
3/20X,«X« ,10X,«YBN« ,//, (12X,2F12o6//) ) 

12    CONTINUE 

DO    11     1=1,1'  1 

YB( I )=YBN( I ) 
11    CONTINUE 

NI=NI+1 

IF    (NI0LTo25)     GO    TO    27D 

GO    TO    <5oo 
2  70    IF    (ERRoGTo (0rrM] ) )    GO    TO    3 

GO    TO    °9Q 
9  99    WRITE     (6,998) 
998    FORMAT     CI'  »2X) 

STOP 

END 


h9 


COMDUTER    PROGRAM    7o 
NEWTON'S    METHOD    APPLIED    TO 
BESSFL.  «S    EQUATION    0<=    ORDER    ZERO 


I  ) 

-.  r 


i ) 


) 

)-(OM**2 
2)*DX/2o 


DIMENSION    X(!" 
l,Y01(lr 1) , YH2( 
DO    1     1=1,1"] 
X(  I)=0    "24rSv(  1-1  ) 

1  CONTINUE 
Y(l)=l    ' 
YP(1  )=0or 
YPN^YP( 1) 
DX=0oP?4    5 
NI=1 

YO(l)=r  „ 
VPQ(1 )=r o 
YPON=YPO(l } 

C    GET    MODF    AND    OMEGA     BY    NE'WT 
C    USE    RUNGA    KUTTA    2ND    OROEP 

OM  =  lo4 
ICO    XN=X(2) 

YN=Y<1 ) 

YON=\C(  1 

YPN=YP(  1 
1 ~o5*0M** 

1  =  1 

YN=Y(1 ) +  (YPP ) +Y 

YPN=YP( 1 )-(OM**2 
1+YPN/XN  )  --DX/Po 

YN=Y(1 )+(YP(3  ) +v 

YPr»N=-YPO(  l)-(  -OM 
l/2o-2o*OM* (Y( 1 )+ 

YON=YO(l )+(vi 0(1 

YP0N=YD0( 1 )-(-0M 
3 ^ YON) *DX/2o-2 , *0 

YCN=Y0(1 )+(YPO(! 

YP( 2)=YPM 

Y(?)=YN 

Y0(2) =YON 

YPn( ? )=Y^n" 

no    2    1=2,11 

J=-l 

YN=Y( I ) 

Y0N=VQ( T ) 

YPN=YP( I )-<0M**2 
DX/2 

YN=Y(I )+(YP(I }+Y 

YPON=YPO( I )-( YPO 

+Y0N)*PX/2o-2o*0 

Y0N=V0(I )+(YPO( I 

1^    (IoGToS)    GO    T 

CONTINUE 

IF     (JoGTc 

XN=X( 1  +  1  ) 

J=l 

GO    TO    ?r 

YP( 1+1 )=vpn 

Y( 1+1  )=YN 

Y0( T+l ) =YON 

ypou  +  i  )=Ypr,r; 

2  CONTINUE 
IF    (NIoLTo*)     GO 
WRITE     (6,1T )    NT, 

ir     FORMAT     ( »1' ,//// 
3  '  ,     OMEGA*     '  ,F12o 
11     CONTINUE 


)  ,  YPir  1)  ,Y0(  iri)  ,YP0(1'--1  ) 


ONS    METHOD*,     FIRST    GUESS    CM    =    2c  P 
INTEGRATION    ROUTINE 


)*(YU  )+YN)*DX/2o-(-o?*(OM**2)*Y(  1) 


pn  )  *nx 

!*(Y(1 

DM )*0X 

-V(  1   )- 

YN  ) 

J+VPON 
-•  Y  (  1 )  + 
M*(Y(1 

)+YPON 


/2o 

)+YN)*DX/2 - -(-o5*(0M**2)*Y( I ) 

/«? 

0 M ) *  0 X / 2o - ( OM* -  2  )  -  ( YO { 1 ) + Y  ON ) *DX 

/2o 

)-  nx/p- 

YPON/XN)*DX/2o-(OM**2)*(YOU) 

)+YN)*DX/2o 

)*DX/2 


1  r-f\ 


)*(Y(I  )+YN)*DX/2o-(YP(I )/X<I)+YPN/XN)* 

YO(I  ) 


54 
300 


4CT 


PN) *nx 
(I  ) /X( 
M*(Y<  I 

)+YPON 

n   54 


/?. 


I) +YPCr /XN)-0X/2o-{ 0M**2)*( 
)+YN)*DX/2o 

K0X/2c 


)    GO    TO    40 


TO     11 
CM, (X( 
////, ' 
5, //,1 


I  ),Y(i)  ,i=i,in,ic ) 

' ,15X,»THIS    IS    PATH    N08 ' , I 

BX,  t x»  ,11  X,  « vi  , //(  ■  <   X,2F12o 


2, 
6//)  ) 
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NI=NI+3 
IF  ADMISSABLE  QUIT 
IF  (ABS( V( I  ]  )  ) 

CORRECT  GMFGA 


LTn (OoOOOOl ) )  GC  TO  15 


IF  NIoGTe25  QUIT 
IF  (NIoLTo25) 

15  CONTIN'JF 
STOP 
END 


GO  TO  100 
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